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Abstract 

The thermodynamic properties such as the specific heat are uniquely deter- 
mined by the second moments of the energy distribution for a given ensemble 
averaging. However for small particle numbers the results depend on the 
ensemble chosen. We calculated the higher moments of the distributions of 
some observables for both the canonical and the microcanonical ensemble of 
the same van der Waals clusters. The differences of the resulting thermody- 
namic observables for the two ensembles are calculated in terms of the higher 
moments. We demonstrate how for increasing particle number these terms 
decrease to vanish for bulk material. 

For the calculation of the specific heat within the microcanonical ensemble 
we give a new method based on an analysis of histograms. 
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I. INTRODUCTION 



Van der Waals clusters have been investigated both with microcanonical and canonical 
simulation methods in the past (see 0-^,0 and references therein). While canonical simu- 
lations are done almost exclusively with the Metropolis algorithm [Q, for microcanonical 
simulations there are the alternatives of doing them by molecular dynamics or with the 
Creutz algorithm |]7,n|. 



The main interest in most publications so far lies on the identification and classification 
of phase transitions in van der Waals clusters. The deficiencies of all simulation methods 
mentioned above occur in the region of the so-called Berry phase. The reason for that is the 
large number of isomers which has to be taken into account for a correct calculation of the 
partition function 0. 

Albeit the fact, that we use for the reason of good comparability Argon clusters as our test 
system, in this work we will not be concerned with the questions mentioned above. Instead 
our main interest here lies in the dependence of the higher moments of the statistical dis- 
tribution functions on the particle number and the differences between microcanonical and 
canonical descriptions. It is well known that for the case of the canonical ensemble both 
the internal energy {E) and the constant volume heat capacity Cy = jAj^{{E — {E))'^) are 
approximately proportional to the number of degrees of freedom (see e.g. [Q]). This imme- 
diately sets up that the relative fluctuations in energy become smaller as the system size 



mcreases. 



Ani? = ^ {{E - {E)Y) / {E) ~ Ar-V2 

Indeed one might see this as the reason for the existence of the Berry phase and one expects 
that its temperature width should decrease with increasing particle number. 

The higher moments of state variables are normally not of much interest in thermodynam- 
ics. The reason is that they can easily be calculated by means of the dissipation-fluctuation 
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theorem from the first moments. For small systems, however, higher moments are a unique 
tool to sensitively study the subtle differences of the same thermodynamic quantity from 
partition functions of different ensembles, since the higher moments explore more sensitively 
the form of the distribution (see [|1],0]). 

The van der Waals clusters explored here are an extremely simple example because only 
one independent variable is to be given, while the volume, the surface and so on adjust 
themselves. 



II. COMPUTATIONAL METHOD 

A system of Ar„ clusters is modelled with the usual Lennard- Jones Potential for the 
interatomic binding, 

with parameters e = lO.SmeV and a = 3.405A. 



A. Canonical Monte Carlo 

For the canonical ensemble the partition function 

N 

^ = Ildxi exp{-(3V{xi, ...,Xn)) (3) 

i=i 

is calculated using the well known Metropolis Monte Carlo algorithm. For optimal per- 
formance we use for the calculation of all thermodynamic quantities the histogram method 
of Ferrenberg et al. ||13| , p!^ ,|5| . For convenience we give a short outline of the procedure. 
If we perform R Metropolis simulations at parameters (3^ and store the Monte Carlo data 
as histograms Np{E) with n^^ being the total number of observations in the i'th run, the 
probability distribution is given through 

PpiE) = Np{E)/np = W{E) expi-^E + (4) 
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where W{E) is the density of states, and F is the Helmholtz free energy. Following Ferren- 
berg et al. the normalized probability distribution can be improved by 

D (E) = ^liN^jE) 
^ Ef=in,exp(-(A-/3)E + /5F,) 

where 

^Dft=exp(/3F,), (6) 

E 

which is simply the partition function. The values of Fi can be determined within an additive 
by selfconsistent iteration of and @. Simply spoken, this method improves ^ 13(E) by 
use of simulation data at other histogram data weighted by the number of observations. We 
slightly improved the method by using fitted probability distributions which are sacrificed 
by tests. 

Now the mean value of any observable 0{E) can be calculated easily as a function of /5. 

(OpiE)) = ^ / 0{E)D,{E) . (7) 

ZjjS JE 

B. Microcanonical Monte Carlo 

The microcanonical partition function 

/n n 
\{dxi W dpj 6[H{xi, ...,Xn,Pi,...,Pn)- E] (8) 
i=l j=l 

is approximated with the microcanonical Monte Carlo algorithm invented by Creutz 
||7| JTT|JT0[| . The algorithm simulates the integral 

f n fE-Vo 

Z= l[dxi dEu 5[V{xu ...,Xn) + Eu-E]. (9) 

J Jo 

where Vq is the minimal potential energy, in our case the potential energy of the best single 
cluster configuration. Ed is an extra degree of freedom called demon, which simulates the 
kinetic energy of the system and is restricted to positive values. As in the Metropolis algo- 
rithm new configurations x' are chosen at random. Here the new configuration is accepted 
if 



(10) 



In this case x' is counted as a new configuration and the demon is set to ^ Ed + AV, 
otherwise the step is rejected. Pictorially the demon might be viewed as a tiny thermometer 
thrown into a large swimming pool. If we denote the cluster system by Ac, the demon 
system by Ad and the combined system by A*^ , the conservation of energy can be written 
as Eq + Ed = E^. The expansion of In Z in a taylor series yields 



lnZc(^°-^D) =lnZc(^°) + 



(11) 



E 



d'lnZc 
dE. 



C j£0 



Eh- 



While the first derivative is easily identified as (3 



dE, 



C J^o 



/3 



(12) 



the higher derivatives are in turn derivatives of (3. The probability distribution for the demon 
energy is now given by 



P{ED) = Cexp{-pED + 



'dE 



El... 



(13) 



In the case that the demon energy is sufficiently small compared to the total energy of the 
system , (3 is indeed the inverse of the demon energy as stated in the hterature . 



{Ed) 



EO-Vo 



dEo EdP{Ed) 



1 



(14) 



C. Multiple normal modes 

In the multiple normal modes (MNM) model, described in detail in [0, we take into 
account several isomers of a cluster, characterizing each isomer by it's binding energy, per- 
mutational degeneracy, and normal modes spectrum. 

The ensemble partition functions are constructed from the single isomer partition functions 
with proper ensemble dependent weights. 
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For the statistical equilibrium of isomers, the calculated one-isomer partition functions 
have to be multiplied by a factor cTj reflecting the permutational degeneracy Rj. In order to 
relate all of them to a common energy (all particles free) also the exponential of the binding 
energy Ei appears as a relative weight between configurations in the canonical partition 
function 

Z = Y:c7,e-^'''Z^ . (15) 

i 

Within the normal modes analysis for a given isomer the potential energy is expanded 
up to second order around the ideal equilibrium position of the isomer. 

V{^) = Y.^ira,)^lT.P'i<^', (16) 

a>l3 ct[3ij 

where be the i'th spatial component of the position of the particle a with respect to 
its ideal equilibrium position, xj^ = 0. Diagonalization of the matrix P yields the 6N-6 
eigenfrequencies [[ 

All investigated isomers were copied from "snapshots" of MC-calculations or constructed 
"by hand" and relaxed numerically until the internal forces were smaller than 10~^'^e/a. 
Then the matrix P was calculated and diagonalized numerically. With ul being the fc-th 
eigenfreqency of the i-th isomer, the canonical partition function is 

^ = E-^e^"' n £7 ■ (17) 
i k=i f^^k 

For the microcanonical partition function at first one has to calculate the micro-canonical 
MNM-partition function for one isomer from the phase space integral: 

Z^E) = I {dp dg)3^-^ 5{e{q,p) - E) (18) 

= 0,N-r.^ n ^ (19) 
fc=i ^k 



^Six eigenvalues vanish due to zero total momentum and total angular momentum. (In the case 
of fully linear structures one has 6N-5 eigenfrequencies.) 
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TABLE I. Total binding energy and relative degeneracies of the most important Ari3 isomers. 



Isomer 



R 



pure icosahedron 



44.33 



1 



singly decorated 



41.47 



180 



doubly dec, neighb 



40.62 



900 



doubly dec, dist. 



39.71 



4800 



where On is the surface of the n-dimensional unit-sphere. To obtain the MNM-partition 
function, all isomer-partition functions have to be related to a common reference-energy, 
and summed up with the relative permutational degeneracies ai , 



The multiple normal modes model gives a unique chance to study the differences be- 
tween micro canonical and canonical ensembles, because for both ensembles all thermody- 
namic quantities can be calculated exactly. Fig.l shows the caloric curves for Ari3 clusters. 
In our model calculations the four most important isomers of Ari3, for which the binding 
energies and relative degeneracies are given in Table 1, have been included. Both curves 
have a similar form but differ significantly by a certain amount of energy. A quite nice fea- 
ture of the model is that the so-called van der Waals loops, often encountered in molecular 
dynamics simulations, can qualitatively be reproduced by insertion of an activation barrier 
into the microcanonical model, i.e. by reducing the accessible phase space up to a specified 
energy barrier. The differences between the ensembles get more drastical if one examines 
the specific heat ( see Fig. 2). The phase transition occurs in both descriptions at the same 
temperature but is much sharper in the microcanonical ensemble. 




(20) 



From ( ]T7| ) and (^Dj) now all thermodynamic properties can be calculated. 



III. RESULTS 
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FIG. 1. Canonical versus microcanonical caloric curve of Aria calculated within the MNM 
model. The inclusion of a virtual activation barrier between the icosahedral and the singly deco- 
rated isomer of about 60 meV in the microcanonical model results in a van der Waals loop which 
is often encountered in molecular dynamics simulations. 
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FIG. 2. Canonical versus microcanonical specific heat of Aris calculated within the multiple 
normal modes model. 
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FIG. 3. Dependence of the relative energy fluctuation on the number of degrees of freedom 
calculated with canonical Monte Carlo at T ~ 5 K. 

Besides the deficiency, that within the multiple normal modes model only harmonic exci- 
tations of each isomer are considered, it is very difficult to find all important isomers and 
their permutational degeneracy as the cluster size increases. In Monte Carlo calculations 
these features are included automatically, but they are exact only within the limit of infinite 
computation time. 

Fig. 3 shows the dependence of the relative energy fiuctuation (1) as a function of the cluster 
size at a temperature of ~ 5 K calculated from canonical Monte Carlo simulations. As ex- 
pected, / ArE is approximately proportional to the inverse of the square root of the particle 
number. 

To get comparable results between microcanonical and canonical simulations a primary 
task is to calculate the temperature within the microcanonical simulations. The simple 
choice to do that is, of course, making use of equation (14) by neglecting the higher order 
terms. Besides that we tried another way, which enables us to calculate the derivatives of f3, 
too by fitting the histogram data of the demon energy to the probability function (|l^). Fig. 
4 displays an example of such a fit. Albeit the difference between the calculated /3's is quite 
small it is not negligible. By chance we have found a cheap way to calculate the derivatives 
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FIG. 4. Fit of the demon energy histogram data to the theoretical probabihty distribution 
function (13). The histogram was obtained for an Aria cluster at E = -445 meV. Prn is obtained 
from the mean value equation (13). 

of f3. 

The caloric curves of our Monte Carlo simulations reveal a behavior very similar to that 
found within the MNM model (see Fig. 5) . Preliminary results show that the relative fluc- 
tuations of the temperature in the microcanonical ensemble decreases similar to the relative 
energy fluctuation shown in Fig. 3 for the canonical ensemble. It is very difficult to com- 
pare the results for the higher moments of the different ensemble quantitatively because the 
differences depend not only on the particle number but also very strong on the temperature 
(see Fig. 2.). Probably other system with a smaller transition phase, e.g. Ising systems, are 
more appropriate for a quantitative study of such dependencies. Nonetheless the decrease of 
the relative fluctuations of the state variables with increasing system size in both ensemble 
indicates that the differences between the ensembles decrease with pa iV~^/^ too. 



IV. CONCLUSIONS & OUTLOOK 



We have shown that in the case of small systems it is not unimportant which thermody- 
namical ensemble is used to describe the clusters. 
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FIG. 5. Caloric curve for Ari3 from microcanonical and canonical Monte Carlo simulations. 
Every data point was evaluated with 8 * 10^ steps. 

Things might get more interesting if more state variables are considered, e.g. by introducing 
spin degrees of freedom along with a coupling to an external magnetic field. A study for 
such a system is in preparation. To investigate the size dependencies more properly not only 
the range of the cluster size (in this study up to 40) has to be expanded. Because of possible 
magic number effects in every size region some neighbor numbers should be considered in 
order to average out such effects. 

The accurate calculation of higher moments, or derivatives of /3 within the microcanonical 
ensemble, remains a problem to be solved. One solution could be to invent a procedure 
similar to the optimized data analysis of Ferrenberg. 
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